{
"cells": [
{
"cell_type": "markdown",
"id": "1b376c36",
"metadata": {},
"source": [
"---\n",
"format:\n",
" html:\n",
" other-links:\n",
" - text: This notebook\n",
" href: L6-7 Solving nonlinear equations in 1d II.ipynb\n",
"---"
]
},
{
"cell_type": "markdown",
"id": "aad6bc09-0de2-4531-9043-0bacc8e41409",
"metadata": {},
"source": [
"# Solving nonlinear equations in 1d II"
]
},
{
"cell_type": "markdown",
"id": "c79162dc",
"metadata": {},
"source": [
"Note. This chapter will be mainly based on Chapter 1 of An Introduction to Numerical Analysis by Suli, Mayers. These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "76076758",
"metadata": {
"jupyter": {
"source_hidden": true
}
},
"outputs": [
{
"data": {
"text/plain": [
"μ (generic function with 1 method)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# | include: false\n",
"\n",
"using Plots\n",
"using LaTeXStrings\n",
"using Polynomials\n",
"using PrettyTables\n",
"\n",
"function simple_iteration( g, x1; N=100, tol=1e-10 )\n",
" x = [ x1 ]\n",
" for n in 2:N\n",
" push!( x, g(x[n-1]) )\n",
" if (abs(g(x[end]) - x[end]) < tol)\n",
" break\n",
" elseif (x[end] == Inf)\n",
" @warn \"simple iteration diverges to Inf\";\n",
" break\n",
" elseif (x[end] == -Inf)\n",
" @warn \"simple iteration diverges to -Inf\";\n",
" break\n",
" end \n",
" end\n",
" return x\n",
"end\n",
"\n",
"function relaxation( f, λ, x1; N=100, tol=1e-10)\n",
" x = [x1]\n",
" r = 0.;\n",
" for n in 2:N\n",
" push!( x, x[n-1] - λ*f(x[n-1]) )\n",
" r = abs(f(x[end]));\n",
" if (r < tol)\n",
" return x\n",
" end\n",
" end\n",
" @warn \"max interations with |f| = $r\";\n",
" return x\n",
"end\n",
"\n",
"function Newton( f, f_prime, x1; N=100, tol=1e-10)\n",
" x = [x1]\n",
" for n in 2:N\n",
" push!( x, x[n-1] - f(x[n-1])/f_prime(x[n-1]) )\n",
" r = abs(f(x[end]));\n",
" if (r < tol)\n",
" return x\n",
" end\n",
" end\n",
" @warn \"max interations |f| = $r\";\n",
" return x\n",
"end\n",
"\n",
"function orderOfConvergence( x, ξ; α=0)\n",
" err = @. abs(x - ξ)\n",
" logerr = @. log10( err )\n",
" ratios = [NaN; [logerr[i+1] / logerr[i] for i in 1:length(logerr)-1]]\n",
" if (α == 0) \n",
" α = ratios[end]\n",
" αr = round(α, sigdigits=3)\n",
" end\n",
" mu = [NaN; [err[i+1] / err[i]^α for i in 1:length(err)-1]]\n",
" pretty_table(\n",
" [1:length(x) err logerr ratios mu];\n",
" column_labels = [\"iteration\", \"absolute error\", \"log error\", \"alpha\", \"mu (α = $α)\" ] \n",
" )\n",
"end\n",
"\n",
"function μ( x, ξ; α=1 )\n",
" return @. abs( x[2:end] - ξ ) / ( abs(x[1:end-1] - ξ )^α );\n",
"end "
]
},
{
"cell_type": "markdown",
"id": "7180faba",
"metadata": {},
"source": [
"## Results from Calculus"
]
},
{
"cell_type": "markdown",
"id": "afcfc0aa",
"metadata": {},
"source": [
"⚠ Intermediate Value Theorem \n",
"\n",
"Suppose that $f : [a,b] \\to \\mathbb R$ is continuous. Then, $f$ is bounded on $[a,b]$ and for all $y \\in [\\min_{[a,b]} f, \\max_{[a,b]} f ]$ there exists $x\\in[a,b]$ such that $f(x) = y$. \n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "7dd70195",
"metadata": {},
"source": [
"⚠ Mean Value Theorem \n",
"\n",
"Suppose that $f \\colon [a,b] \\to \\mathbb R$ is continuous on $[a,b]$ and differentiable on $(a,b)$. Then, there exists $c\\in[a,b]$ for which \n",
"\n",
"\\begin{align*}\n",
" f(b) - f(a) = f'(c) (b-a)\n",
"\\end{align*}\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "352888f3",
"metadata": {},
"source": [
"*Proof:* Consider $g(x) := f(x) - \\frac{f(b)-f(a)}{b-a} x$, a function for which $g(a) = g(b)$. You can show that there exists $c\\in [a,b]$ for which $g'(c) = 0$ (this is known as Rolle's theorem). This follows from the intermediate value theorem - either the max or min is achieved at $a$ and $b$ (and the function is constant), or there is an extreme point in $(a,b)$. In the latter case, you can show this point is stationary: i.e. $g'$ vanishes here. \n",
"\n",
"Cauchy's Mean Value Theorem (need this for the proof of Taylor remainder theorem)\n",
"\n",
"Suppose that $f, g \\colon [a,b] \\to \\mathbb R$ are continuous on $[a,b]$ and differentiable on $(a,b)$. Then, there exists $c\\in[a,b]$ for which \n",
"\n",
"\\begin{align*}\n",
" \\big( f(b) - f(a) \\big) g'(c) = \\big(g(b) - g(a) \\big) f'(c)\n",
"\\end{align*}\n",
"\n",
"*Proof:* $h(x) = \\big( g(b) - g(a) \\big) f(x) - \\big( f(b) - f(a) \\big) g(x)$ is such that $h(a) = h(b)$ and so there exists $c \\in [a,b]$ such that $h'(c) = 0$."
]
},
{
"cell_type": "markdown",
"id": "50a9ac6d",
"metadata": {},
"source": [
"⚠ Taylor Remainder Theorem \n",
"\n",
"Suppose that $f\\colon [a,b] \\to \\mathbb R$ is a function with $n$ times continuously differentiable. Then, there exists $c\\in[a,b]$ for which \n",
"\n",
"\\begin{align*}\n",
" f(x) = f(a) + f'(a) (x - a) + \\dots + \\frac{f^{(n-1)}(a)}{(n-1)!} (x - a)^{n-1} + \\frac{f^{(n)}(c)}{n!} (x-a)^{n}\n",
"\\end{align*}\n",
"\n",
"
\n"
]
},
{
"cell_type": "markdown",
"id": "7a421ef6",
"metadata": {},
"source": [
"\n",
"*Proof:* Define $F(x) = f(x) - \\sum_{k=0}^{n-1} \\frac{f^{(k)}(a)}{k!} (x - a)^k$ and $G(x) = (x-a)^n$ and note that $F(a) = F'(a) = \\cdots = F^{(n-1)}(a) = 0$ and $G(a) = \\cdots = G^{(n-1)}(a) = 0$ and $G^{(n)}(x) = n!$. Applying Cauchy's Mean Value Theorem $n$ times, we can conclude there exists $c_1,\\dots,c_n$ between $a$ and $x$ such that\n",
"\n",
"\\begin{align*}\n",
" \\frac{F(x)}{G(x)} &= \\frac{F(x) - F(a)}{G(x) - G(a)}\n",
" = \\frac{F'(c_1)}{G'(c_1)} \\nonumber\\\\ \n",
" &= \\frac{F'(c_1) - F'(a)}{G'(c_1) - G'(a)}\n",
" = \\frac{F''(c_2)}{G''(c_2)} \\nonumber\\\\ \n",
" &= \\cdots = \\frac{F^{(n)}(c_n)}{G^{(n)}(c_n)}\n",
" %\n",
" = \\frac\n",
" {f^{(n)}(c_n)}\n",
" { n! }\n",
"\\end{align*}\n"
]
},
{
"cell_type": "markdown",
"id": "08a20cbb",
"metadata": {},
"source": [
"## Previously...."
]
},
{
"cell_type": "markdown",
"id": "a0fc8616",
"metadata": {},
"source": [
"⚠ Change of Sign Theorem \n",
"\n",
"Suppose that $f\\colon [a,b] \\to \\mathbb R$ is a continuous function with $f(a) f(b) \\leq 0$. Then, $f$ has a *root* (or *zero*) $\\xi \\in [a,b]$ (i.e. $f(\\xi) = 0$). \n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "e7e37d68",
"metadata": {},
"source": [
"*Proof.* If $f(a) f(b) = 0$ then $a$ or $b$ is a root of $f$. Otherwise, $0 \\in (\\min_{[a,b]}f, \\max_{[a,b]}f)$ and so the Intermediate Value Theorem applies. "
]
},
{
"cell_type": "markdown",
"id": "61ef8668",
"metadata": {},
"source": [
"⚠ Brouwer's Fixed Point Theorem. \n",
"\n",
"Suppose that $g \\colon [a,b] \\to [a,b]$ is continuous. Then, there exists a fixed point $\\xi \\in [a,b]$.\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "0c965102",
"metadata": {},
"source": [
"*Proof.* We can apply the change of sign theorem to $f(x) := g(x) - x$."
]
},
{
"cell_type": "markdown",
"id": "c62bc0cc",
"metadata": {},
"source": [
"Definition. (Contraction) \n",
"\n",
"A function $g : [a,b] \\to \\mathbb R$ is a *contraction* if there exists $L \\in (0,1)$ such that \n",
"\n",
"\\begin{align}\n",
" |g(x)-g(y)|\\leq L|x-y|\n",
"\\end{align}\n",
"\n",
"for all $x,y \\in [a,b]$,\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "766b8b5f",
"metadata": {},
"source": [
"By the Mean Value Theorem, if $|g'|\\leq L < 1$ on $[a,b]$, then $g$ is a contraction on $[a,b]$. "
]
},
{
"cell_type": "markdown",
"id": "32f25a86",
"metadata": {},
"source": [
"⚠ Contraction Mapping Theorem (also called Banach's fixed point theorem) \n",
"\n",
"Suppose $g:[a,b] \\to[a,b]$ is a contraction. Then, there exists a unique fixed point $\\xi = g(\\xi) \\in [a,b]$. Moreover, the iteration $x_{n+1} = g(x_n)$ converges at least linearly to $\\xi$ for all $x_1 \\in [a,b]$.\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "e240a44e",
"metadata": {},
"source": [
"*Proof.* Existence of a fixed point $\\xi = g(\\xi)$ follows from Brouwer's fixed point theorem. If there exists $\\zeta = g(\\zeta) \\in [a,b]$ then $|\\xi - \\zeta| = |g(\\xi) - g(\\zeta)| \\leq L |\\xi - \\zeta|$. Since $L \\in (0,1)$, we must have $\\xi = \\zeta$ and the fixed point is unique. \n",
"\n",
"Notice that $|x_{n+1} - \\xi| = |g(x_n) - g(\\xi)| \\leq L |x_n - \\xi|$ for all $n \\geq 1$. Therefore, we have $|x_{n+1} - \\xi| \\leq L^n |x_1 - \\xi|$ which goes to zero as $n\\to\\infty$ (as $|L|<1$). Finally, we note that by the Mean Value Theorem there exist $\\eta_n$ between $x_n$ and $\\xi$ such that\n",
"\n",
"\\begin{align}\n",
" \\frac\n",
" {|x_{n+1} - \\xi|}\n",
" {|x_n - \\xi|}\n",
" %\n",
" = \\frac\n",
" {|g(x_{n}) - g(\\xi)|}\n",
" {|x_n - \\xi|}\n",
" = |g'(\\eta_n)|.\n",
"\\end{align}\n",
"\n",
"This quantity converges to $|g'(\\xi)|$ as $n \\to \\infty$.\n"
]
},
{
"cell_type": "markdown",
"id": "3623a6c4",
"metadata": {},
"source": [
"Theorem. \n",
"\n",
"Suppose $g$ is continuously differentiable with fixed point $\\xi$. Then, if $|g'(\\xi)| >1$, then $\\xi$ is an unstable fixed point of $g$.\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "63a6e0c9",
"metadata": {},
"source": [
"*Proof.* Since $g'$ is continuous, there exists an interval $I$ containing $\\xi$ such that $|g'| \\geq L > 1$ on $I$. Consider the iteration $x_{n+1} = g(x_n)$. By the intermediate value theorem, if $x_1 \\in I \\setminus\\{\\xi\\}$, there exists $\\eta_1 \\in I$ for which $|x_{2} - \\xi| = |g(x_1) - g(\\xi)| = |g'(\\eta_1) (x_1 - \\xi)| \\geq L | x_1 - \\xi |$. Similarly, if $x_{2}, \\dots, x_{N} \\in I$ then \n",
"\n",
"\\begin{align}\n",
" |x_{N+1} - \\xi| = |g(x_{N}) - g(\\xi)|\n",
" %\n",
" &\\geq L|x_{N}-\\xi| \\geq \\cdots \\geq L^{N} |x_{1}-\\xi|.\n",
"\\end{align}\n",
"\n",
"Since $L>1$ and $x_1 \\not=\\xi$, we have $x_{N+1} \\not\\in I$ for sufficiently large $N$. If the sequence returns to $I$: if there exists $m$ such that $x_m \\in I$, then by the exact same argument as above $x_{m+M} \\not\\in I$ for sufficiently large $M$. As a result, $(x_n) \\not\\to \\xi$. "
]
},
{
"cell_type": "markdown",
"id": "e61f4cfb",
"metadata": {},
"source": [
"## Relaxation\n",
"\n",
"* Define $x_{n+1} = x_{n} - \\lambda f(x_{n})$ "
]
},
{
"cell_type": "markdown",
"id": "fd1cca2e",
"metadata": {},
"source": [
"Example (Kepler's Equation) \n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "966f7a90",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dd3gU1d4H8N+UZNN7L/QaepcOIQmEIog0FZCr2PXaEAs2wKu+V5R7vejFio1LVRCQFlroIDUQqgRDGqmQvpudmfP+sZsQU0jb7Gz5fh4fn9nZk52TYXa+OWfOnOEYYwQAAGCveLUrAAAAoCYEIQAA2DUEIQAA2DUEIQAA2DUEIQAA2DUEIQAA2DUEIQAA2DUEIQAA2DUEIQAA2DUEIQAA2DXLCkLG2Msvv1z/8nq9vvkqA3enKIosy2rXwn7h4FcRdr6KZFk2+cygnEXNNSpJkrOzc/0PsqKiIjc3t2atEtRGr9criqLRaNSuiJ0qLCx0d3dXuxZ2CjtfRVqtVhRFURRN+JmW1SIEAAAwMwQhAADYNQQhAADYNQQhAADYNQQhAADYNQQhAADYNQQhAADYNQQhAABYB0mhb69y6/808e3vprwnEQAAoJkwotnx8qprPBF5ObExYZypPhktQgAAsALzj8mrrimGZd5kIUhkGy1CRVEeeeSRkpIStStiblOmTJk2bZratQAAaHb/Pq8sOWdMwac6UUyoKZPQFuYa1el07u7uK1eubNa6WZqdO3c6Ojp+9tlnalUAc42qC9Ndqgg738zWX1em75EVRkQ0MZytieQ1DqZsxdlCi5CIeJ6fOnWq2rUwq+zs7MTERLVrAQDQvI5ksVn7jCk4OJBbMVgWOBNf1MM1QgAAsFDXC9mkOEkrExF18uJ+jRGdm6H5hiAEAABLdLuMxu+Qs0qJiAKcaetowbd5LsUgCAEAwOLoFZq2W7pwmxGRk0AbosTW7iYdKloJghAAACzO80fkuDRGRBzR10OFQYHNlYKEIAQAAEvzUYLy34vGmyUW9xUeate8UWUjo0Yt2W+//Xb79u2+fft27NhR7boAAFi631LY67/LhuUZbfk3ejZ7gw0twuaVm5v7448/pqWl5eXlqV0XAABLdyqHTd8tyYyIaGgQ990woRm7RMuhRdi8Dh482K1bt/nz56tdEQAAS5dewibGycUSEVEbd+7nKFEjmGO7aBE2oz179qxYsSIlJeX48eNq1wUAwKIV6mnsdjm1mBGRj4a2jRH8ncy0adtsEW74U/nikiKba/I4R56e7ypUn/suMjLyk08+efvtt0NCQsxUFQAAKyQzemivfDaPEZEDT+ujxA6eZugTNbLBIDQ8qqOovvOVmsb5W3LyjBp2ZlpaGlIQAODuXjgib76hUPnNEiODzZeCZJNdoxxRdx+z7kQi6lHTFnNycgICAsxcEwAA6/JporLsgvFmibd68bPbmzuYbLBFSES7YsWj2UxWzLQ5R4HuCaghCE+fPt2zZ08zVQIAwAptS2EvHzPeLDG1Nf9uH7MMj/kr2wxCZ5HM3LKu0ZkzZ+655x61awEAYKESb7EH9kqSQkTUz5/7brg5bpaozga7Ri1BVlbWH3/8ceHChVGjRqldFwAAS5RRQrHb5fwyIqJW7tyWGNFFpaaZbbYIVffTTz/p9fqePXt6eHioXRcAAItTKtGkOCmlmBGRhwNtjhECnFWrDIKwWbz44ospKSktWrRQuyIAABZHYfTgXvl4tvFmiZ+jxK7eal7MQtdos+A4DikIAFCjV47LG5ONoxk/HShEVbsJ28zQIrRcZ8+eTUlJ4Xk+NjaW49Qf+wMA0HRfX1Y+OWdMwfnd+Sc7q98eU78GUKOEhIR169aNHz8+ODh4yZIlalcHAMAE4jPYM4eMN0tMbsV/0E+FmyWqQxBaqK+//jomJoaIevXq9b///U9RzHVTJABA80gqZFN2S2UKEVEfP+7HEQJvGV1dCMJmIUnS5s2bFy9ebHhZsVB/p06dqhhxWlpamp6ebsr6AQCYV6GeJsXJOVoioiBn2hAtqHWzRHUWUxGTKjm5t+jgZpJMMN8o5+DoNuJ+5+6DGvRTW7ZsiYyMfOedd958882ioqK4uLi33nrL8FZmZub+/fur/8jAgQPDwsIqXubn52s0GsOyi4tLXl5e5XcBAKyIwuihvfK5PEZETgJtiBbDXS2jMUhEthmEjN1a8y9WpjPV50nZ6Q0NwtjY2PPnz/fu3ZvjuEOHDg0ZMqTircDAwKlTp9b5Ce7u7qWlpYbl4uJid3f3BlUAAMByvHr8L3Nq1zgnpYpsMQg5zrFlJ93Vs6b6PMdWnRr6IxqNZuvWrRMmTCCiffv2RUZGVryVnZ29b9++6j9SpUXYrl27goICw7Isy3iEBQBYqR+uKkvKh4m+0ZN/qJ3FXZKzxSAk8nvyH1JGMmMmGGDC8bxDcKtG/GBJSUloaGhpaelvv/1W0S9KRP7+/vVpEc6aNevkyZMjRoxISUkZOnRoRTcpAIAVOZzJHj9oHCY6qSW/SI05tetkm0HICaJDWFt16zBv3rzt27cnJycHBAS4uro29Mejo6O1Wu2uXbuSk5OXLVvWHDUEAGhWyUVs8i5JJxMR9fCxoGGiVdhmEKru8uXLn3zyyRdffPHPf/7zlVdeadyHGHpWAQCsUZGeJu6UM0uJiPyc6Jdowc1B7TrVwuL6am2Dh4fH4MGD169f37FjxzFjxqhdHQAAs2JEjx6Qz+YxInLkaf0osY27RTYGiQgtwmYSHBw8e/ZstWsBAKCOBb/La5OMozQ+GywMt4AHxN4FWoQAAGBK664rH569M5vo3I6WHjSWXj8AALAiJ3PYnHiZERHRmDDufcuYTfTuEIQAAGAa6SVsYpxcIhERdfbiVkeKgkX3iRohCAEAwARKJZoUJ6cVMyLy1dDmGMHTUe061Q+CEAAAmooRPXJA/r38ofProsS2HtbQGCQiBCEAADTduyfl1deMA2T+M0gYadnDRKuwhdsneJ6XZblv375qV8SssrOzJ02apHYtAADo5+vK4tPGFHyxK/9EJytrYtlCEDo4OCQmJhYWFqpdEXPr0KGD2lUAAHt3MofNrjRM9KMBVjBMtApbCEJCJAAAqMFKh4lWYWUNWAAAsBAlEk3caRwm6qOhTdYzTLQKBCEAADSYYTbREznlw0RHie2sZ5hoFQhCAABosHcqDRP9dKAQGWKtKUgIQgAAaKh115X3Kg0TfbKzdUeJddceAADMrPJsoqOtc5hoFQhCAACoryrDRNdY5zDRKhCEAABQL6WVhon6WvMw0SpMcx+hVqu9ePGioiidO3d2cXGpXiAjI6O0tNSw7ODgEB4ebpLtAgCAeSiMZuw1DhPVCLQx2oqHiVZhgiDctGnTnDlzWrVqJQhCSkrK6tWrR4wYUaXMzJkzExMTXV1diahNmzZxcXFN3y4AAJjNK8flTcnGATJfDBGGBNlICpJJgrB9+/bnzp0LDQ0log8++ODZZ589f/589WLLli2bMmVK0zcHAABmtvyi8sk5Ywq+2oN/uL1NXVYzwS/TuXNnQwoSUd++fbOzs2ssVlBQ8Mcff+j1+qZvEQAAzGZnGnvuiGxYntyKf7+v1Q8TrcLEc41+9tlntTX7Fi5cqNFobt68+cEHHzzzzDO1fQJjrHLHacuWLdu1a1dbYUVRFEVpSoWh0ZRyalfETmHnq8iudv7F22z6bkVSiIj6+HHfDeOIKQpTrT4NPfPwfN3tPY4xk/1CixYtWrdu3aFDhzw8PKq8lZOT4+fnR0QHDx6MiYk5cOBAnz59qn+CJEmOjo6VLzGOHTv2ySefrG2LxcXFhuuOYH56vV5RFI1Go3ZF7FRRUZGbm5vatbBT9rPzc3XcyDiH60UcEYW40N4oXUgNoyHNSqvViqIoivVtxTk5OdVZ2GQtwiVLlqxcuTI+Pr56ChKRIQWJaMiQIQMHDjx06FCNQUhEgiDs2bOn/tu1k8PRAiEI1cUYw8GvFjvZ+VqZHtgjXS9iROTuQNvGiB18HNSuFInlTPmZJvmUTz/9dNmyZfHx8UFBQXcvKctySkpKRS4CAIAFYkSP7pePZDEiEjhaOVLo7mM7w0SrMEEQbty48YUXXnjqqadWr15tWDNv3jxBEN57773ExMRVq1ZlZWUtWLBg2LBhGo3mp59+UhRlwoQJTd8uAAA0kzdPyP8rn1P73wOFCS1saphoFSYIQm9v7/nz5xPRrVu3DGsM1x379+9vuHHe3d29devWO3bsIKKBAwd+//337u7uTd8uAAA0hy8vKe+fMabg37vwz0TYcgqSaQfLNJ0kSc7OzvW/xcJ+LllbIFwjVFdhYSH+oFSLbe/831LYpDjJMEx0fAtuY7RlzSba0MEy9WHjOQ8AAPV3MofN2CNV3CyxaqRlpWAzQRACAAAR0fVCNn6HVKQnImrtzv02WnRTf5SoOSAIAQCA8nQ0dod8s5SIyFdD28YIgc5q18lcEIQAAPZOK9PEOOnSbUZETgJtjBY7etpBl2g5BCEAgF1TGM3aJx+8yYiI5+inETb1ZIn6QBACANi1l4/J668bb5b4eIBwf2u7ywW7+4UBAKDCfy8q/zpvTMEXu/IvdLXHULDH3xkAAIjo5+vKs4eNz1ea0ppfMsDWnq9UTwhCAAB7tCedzdwnGx6oNDiQ+3GEwNvXlcE7EIQAAHbn92w2KU7SykREHT25X2NEJzttDRIhCAEA7M3VfDZhp1SoJyIKdeW2jxF87XuqRAQhAIAdSS1m0dvkzFIiIj8n2hkrtHK31y7RcghCAAB7kaOl6G1ychEjIheRfo0WI7zsPQUJQQgAYCcK9DRmu3H6GEeefokSBwUiBYkQhAAA9qBUont3SidzjE+c/2mkMDoMKWiEIAQAsHEyo5n75PgMRkQc0fIhwlT7mz7mLrAvAABsGSN6/ID8y5/G6WP+r78wtyPO/H+B3QEAYMteOCJ/e8WYgq/14F/pjtN+VdgjAAA2a94x+dNEYwo+3ol/v58d3zZfOwQhAIBtevW4/PE5YwpOb8N/PljA8JgaIQgBAGzQmyfkfyYYU3ByK/7HEQJisDYIQgAAW/P2SfkfZ4wpeF8rfnWk4ICTfe2wbwAAbMo7J+XFp40pGBvOrRqJFKwDdg8AgO1495S8qDwFx4RxG6JEDcbH1AVBCABgIz5KUBaeMqbg6DBuQzRSsF4QhAAAtmDJOWX+cePj5mNCuY3Rdv2IwQZBEAIAWL3Fp5VXjhlTcHSYvT9ot6FEtSsAAACNx4hePiovPW/sEY0O5TZEIQUbBkEIAGCtZEaPHZBXlM+gZugRdcZ5vYGwwwAArFKZQjP3yuuuG1NwYkt+daSAtmAjIAgBAKxPsUT375J2pDLDy1nt+G+HCSJGfTQKghAAwMrcLqNxO6TDmcYUfK4L/++BmECt8RCEAADWJLOURm+TzuYZU/DVHvyHeKZE0yAIAQCsxp+FLGa7fDXf+Kz5T+4RXuiK/tCmQhACAFiH07ls/A45vYQRkQNPK4YJD7VDCpoAghAAwAr8mqw8tFculoiInEVaGymOb4HLgqaBIAQAsHRLzimvHpcVRkTkraENUeLwYKSgySAIAQAsl6TQ80flzy8YbxZs485tGS109kIKmhKCEADAQt3S0ZTd0p504wDRQYHchigxwFndStkgBCEAgCW6VsAm7JQv3jam4Iy2/IphmDimWSAIAQAszqFMdl+clK0lIuKI3u7Nv9sbGdhcEIQAAJbly0vKc4flMoWIyEWkH4YL97fGbRLNCEEIAGApivT01CH5pz+MQ2OCXWhTjNjXD0NjmheCEADAIiTeYlN337ko2NOX2xQjhLsiBZsdghAAQH0/XFWePmS8X56IZrXjlw8RXHCGNgvsZgAANWllevW4/GmisTvUWaRPBwpzO+KioPkgCAEAVHM5n03bLSeUP0qikxe3NlLo5oPuULNCEAIAqGPlH8qTh+QivfHl7Pb854MFV5yVzQ67HADA3LJKad4x+cfy0aFOAn3YT3geD1RSCYIQAMCsVl9Tnjsi52iNLzt6cutGoTtUTQhCAAAzuVlKTx+SN/xpbAhyRI925JfeI7g5qFsve4cgBAAwh3XXlacP3WkItnTjvhoqRIeiIag+BCEAQPOq3hB8rBO/ZIDgjoagZUAQAgA0F71Cyy8qb5+Ub5cZ17Rx574eJozEY3UtCYIQAKBZbE1hLx+TL5VPmcZz9EwE/0E/3CBhcfAPAgBgYom32MvH5B2prGJNe0/um6HC0CA0BC0RghAAwGTydLTwlPz5RUUyXhAkNwd6uRv/Wg88U9dyIQgBAEygTKH/JCrvnb5zOVDg6LFO/KI+gr+TqjWDuiAIAQCaRK/QqmvKotPKtYI7faGRIdzSe4TuuE3eGiAIAQAaqUyhby8rH5xVbhTdicCOntySAcL4FohAq4EgBABosAI9fXVJ+fd5JaX4TgQGONOCnsJTnXkHTBpqVRCEAAANkFzElpwRf0jSF+jvrAx0ple6C0915vEoXWuEfzQAgLopjHalsy8uKptuKJJy58wZ5EyvdBeeRARaM/zTAQDcTWYpfXdF+eryX8bCEFEXb+6lbvxDbXkN7ouwcghCAIAalEi0MVlZ+YeyM41V3BRIRBzRiCBlfk/H0WEcxsPYBgQhAMAdxRLFpSm/XGcbkpUi/V/e8tHQw+35xzvxoUKRu7tGpQqC6SEIAQAorZhtSWGbkpU96Uwr/+UtjmhIEPdYJ35qa94wO0xhoSp1hOaCIAQAO1Wm0IlsFpfGNt9QTuUwVq1AhBf3YDv+obZcK3d0gtoyBCEA2JESiY5msf03lf0Z7Fg2K5FqKNPNhxsfzk1tw/fyRf7ZBQQhANiyEokSb7GzeSwhj53IZidymF6poZgjT8OCuQkt+AktuNZo/9kZ0wShJEn79u0rKCgYOnSov79/jWXOnTt36dKlzp07d+3a1SQbBQBrxBjbsnnzlvVrLl+6FBQYeM/wkY88/oSHh4epPv96ITuXxxLyKCGPnc1j1wqYXL3Ts1x7T25oIDc6jBsdxns6mqoKYGU4VkPHeMNIkjR69OiCgoK2bdvu3r17586dvXr1qlJmyZIln3zySUxMzM6dO1955ZUXX3yxto9ydnbW6/U1vltdUVGRm5tbk2oPjaXX6xVF0Wgwdk4dhYWF7u7uateiwcrKyqZPnOCee2NaG692Pq7ZJWXxqfnrbhSu2rglIiKiQR91S0dJhSypkKUXU0YpSyqgpEJ2JZ8V1nX+aOPORYVygwO54cFcS7fGNP6sdOfbBq1WK4qiKJqyO9MEQfjzzz+/+eabZ86c0Wg0ixcv/v333zdt2lS5QH5+fmho6OHDh7t3756QkDB48OC0tLQa/wBEEFoRBKG6rPRc/NqLzzsl7JvbLbjyyj9uFb90PPPo2fOOjsZGWYGesktZjpZydJSjZblaytGyLC3lailHx3K0lFzESmu6vFedyFN7D667D9fDl+vuww0M4HyafMxa6c63Dc0RhCb4rF9//XXSpEmGE+L06dMXLlyo1+sdHBwqCuzatatFixbdu3cnou7du4eFhe3Zs2fSpEk1flr3AI/iw1vruekyna64thOxIDp3vYd3rbu/henLShMOMV1pPTd6Nxynad9T9AuuuyRjpYnHlII8E2yUyCG0jWPLTvUpWXb9gj7jT5NslLl5OXTsU5+SUnaa7upZk2yUc3Jx7j6YEx3qLKkU3S5NPEayXGfJuuFYqkkjjqXSMv1v69f8NrlHlfXtvF17ezn2XbSZek80xF5ZTZfxiKhNWfrAooSWRHc58lxFCnfjQl0o1JULc6UQF844BXYxUTFRChUTUdOOpTKtttipUc8YxLFUkwYdS3q9XhAEnq95XnPBy88pon89P6qCCYIwLS1twIABhuXw8HBZlm/evBkeHl65QOWX4eHhqamptX1aVCvfW2s/rf/WS2p/qyisnd+L/67zE26v+Xfpid313+Ld8c5u/m9+yzu53r1Y8YFNBRu/MNVGieP8nl/qEN7+7qXKrl/I/Ww+NbkPoILrpCfFoRPuXkYpLsha8qxpvs9EROQyMNZzyrN1Fsv+z3wp84apNmqBx5Isy7Isk1UdS1dyizq7C3xN87EMD3Rae/pQYat77/LjPnLB1j9ecFUadiwV1f5WU46lRh/QFngsVbCiY+kuPKc86zIwtuIlz/NcXVMAmSAIDflsWDYslJWV1VaAiERRrH/nZ1PIhbfrsyG58JYJN6roSvQlxbxQx2V3Kd80f3MZMVaWn0dBdfyy+oI8Ex5tRKQvuFXnHlaKi1iZ1oQblfLz6vXPWpRvwo1a4LGk1+sNVTLnsVSop+tFXFIh6S/eGtnwY6lMUUSh5j/kRZ7jpDvHiatIvhryd2J+TuTjSL4a5qvhWuhLXS7hWGoMez4vOTg4VA6gGpkgCIODg7Ozsw3LWVlZRBQSElJbAUOZ4OBaG+nnsotcB42t56ar9MFWxvG8y4AYx3p0X4jTny+K38D0ZXWWrA+nzv2c/evugtCMeVAURbnYNF8wh+BWbt0HUl1/9Tj1GsYX5ukzU0yyUc7V02nYJKc693BwODdngfbyKdNs1NHJfcR9Qj3+Wf0fW1jy+y7Gaulia9BGLfJY0uv1hp3ffMdSajE7ncvO5tKVfHatkF0rYJnl7SCe7nkkeG47Xa1dOzXSuWjPxS+t8a3EW9onJvR5aLLoqyE/J86p5hNXWKmrRRxLdznz1LFRizyWKljFeUmWZY7jau0adfd2HzWVc2xYx7UJBsssW7Zs3bp18fHxRPTdd999/vnnx48fJyKdTicIgiiKycnJHTt2TE1N9fPzy87ODg8Pv3r1auXO0goYLGNFMFhGXSYfr8GI/shnp3PZ6Vx2KoedzmXZDWl9BThToDMX7kp+TpyPhnw0hv/fWfB14rwcafrE8eOFzJEtfCv/bIFOmrHzSvyJs97e3ib8jZoPBsuoyEJHjRYUFHTt2nXcuHFdunRZtGjR8uXLJ0+eTESDBw++77775s2bR0SzZ89OSkqaPXv2Dz/80K5du++++67Gj0IQWhEEobpMci4u0lP8TbYrTTmVw87ksoK6vnkagdq4c209qJ0H18adC3WlIGcuzJWCXDjH+j2TPTMzc+zIYQ+1dJncPkDkOSI6l1Xwzu9pCz7618RJ9zXx1zEbBKGKLDQIiSgjI+Obb77Jzc2dOHHiiBEjDCvXr1/ftm1bwz2FkiR9++23iYmJ3bp1mzNnTm2/A4LQiiAI1dXoc7HC6FQu25nK4tKUw5mstvGZROTpSD19uV6+XFdvrq0H19aDwlxN8OChgoKC9995e8/O7ZJOS4LQvmPHBe99aBhVbi0QhCqy3CA0FQShFUEQqquh5+KbpXTgprIrjW2+oWTUMtjay5G6eHN9/Iz/dfbieMw1VhMEoYos9D5CALBYl26zVdeUddfZxds1/MnLEfX05aJDuaFBfC9fCnVF7oE9QhAC2KDUYrYmif3vmnIqp4b8C3HhokO5mDAuKoQPcDZ/7QAsC4IQwHYU6GnNNeV/15T9N5ny1wR0EWlYEBcdyseEcV290fIDuANBCGALLt5myy4oP15Vqkw57STQ+Bb8g2252HC+lpvzAOwdghDAismMNt9QliUqe9L/MuxN4GhUCPdgO/6+VrxHY+78BrAjCEIAq3SrjPv8rPLfi0py0V/6QLt4c4934qe34QNx8Q+gfhCEAFYmR0vLLsj/OueYr7/zPASeo7Hh3PNdhFGhJrjVD8CuIAgBrEZ6CfsoQfnyklIiEZEx7/ydaG5H/qkIPhw3PwA0CoIQwArk6uj9M/LnFxRtpQcsdvPhXu7GT2+DUTAATYIgBLBoxRL967zyUYKcX+k5BH38uJc76WZ0QhsQwAQQhAAWSmb09WVl4Sm58oxo9wRwb/USxoZzhYWlSEEAk0AQAliigzfZc0fkM7l3RoRGeHH/6MdPalm/pzwAQL0hCAEsS0YJvXpc/umPOzPDhLlyb/XiH+nAiwhBgGaAIASwFDKjf51X3j0lF5XPDuMq0us9hZe68s74pgI0G3y9ACzCmVz22AH5RKU5sqe14ZcMwE0RAM0OQQigslKJFp2Wl5xTpPJn5Hbz4T4dKIwIRgQCmAOCEEBNR7LYw/Hy1XxjQ9BZpHd6CS93w+VAAPNBEAKoQ1LovTPyP87caQiODOa+HCq080BDEMCsEIQAKkgqZA/tlY9mGRuC3hr6qL/wSEceGQhgfghCAHNbd115/IB8u3ymmJHB3PcjBAyKAVALghDAfPLL6JnD8so/jJ2hIk8LevJv9RIEhCCAehCEAGayL4PN2ienFhu7Qzt5cStHCL39kIEAKkMQAjQ7RvTxOeX132XDuBiO6InO/McDBBd8/wAsAL6IAM2rSE+PHpDXJhm7Q/2d6Ouhwr2YMhTAYiAIAZrR1Xw2eZd8/paxO3RIELc2Ugx2UbdSAPAXCEKA5rI1hc3cJ93SGV8+3on/zyDBEU1BAAuDIAQwPUb0z7PKGydkwyMknAT6bLDwSAdkIIAlQhACmFh+Gc3cJ225YewObeXO/RIl9PLF6FAAC4UgBDClPwvZ+J1yYvlFwehQblWk6KtRt1IAcDcIQgCTOZbFJsZJmaVERBzR/B78P/riZnkAS4cgBDCN9deV2fFyqUREpBHo22HCg21xURDACiAIAUzg/84qr/8uG/pD/ZxoQ5Q4JAgtQQDrgCAEaBKF0fNH5GUXjPfLd/TktozGo5QArAmCEKDxdDLN2ievu25MwRHB3M9Rog+GxgBYFQQhQCMV6GnSTmlvhnGA6Iy2/PfDcb88gPVBEAI0xs1SGrtdOp1rTMFnI/h/DxTwXF0Aa4QgBGiwPwtZ1Db5WgEjIo7o/X7Caz3QEgSwVghCgIa5nM+ithofKyjy9NUQYQ7mTgOwZghCgAY4m8ditklZpURETgKtGyWOb4H+UADrhiAEqK8TOWzMNilXR0TkKtKGaDE6FCkIYPUQhAD1Ep/BJuyUCvVERN4a2jpavCcAKQhgCxCEAHXbm8Em7JCKJSKiAGfaMUbsiadJANgKBCFAHXans3t3SiUSEVGoK7crVujkhRQEsB0IQhzUzk4AAB3kSURBVIC7ic9gE8tTMNyV2zMO06cB2BoM+wao1c40FlveIxruyu1FCgLYIrQIAWq2K41N3ClpZSKilm7cnnFCG3ekIIANQhAC1ODgTTYpzpiCrdy5vWOFVkhBABuFrlGAqk7ksHHlPaIt3JCCADYOQQjwF+fy2JhtUoGeyHinBFIQwMYhCAHuuJrPRm83zh3j50R7xoq4UwLA5iEIAYxuFLHobXJGCRGRpyNtGyN28UYKAtg+BCEAEVFaMRv5m5xcxIjIRaTNMWJfP6QggF1AEAJQtpZitslJhYyInEXaMlocGoQUBLAXCEKwd/llNGa7dOE2IyIHntZGiiODkYIAdgRBCHatVKJxO6RTOYyIBI5+HCHg+YIA9gZBCPZLUmjaHulQJiMinqMVw4XpbfCNALA7+NqDnWJEjx2Ut9xghpf/ukeY1Q5fBwB7hG8+2KnXf5e/u6IYlt/sxT/XBd8FADuFLz/Yo88vKP931piCs9rxi/oI6tYHAFSEIAS7s/qa8twR2bA8vgX37TABw2MA7BmCEOzLnnQ2Z7+sMCKiAQHc6khRxJcAwL7hHAB25GQOmxQn6WQioggvbuto0RUPIgOwewhCsBfXCti4HVKhnogo1JXbNkbw0ahdJwCwAAhCsAvZWhq9Xc4sJSLy1dDOWKGFG64MAgARghDsQalE9+6UrhWUT6g9WozAw5UAoByCEGycwmhWvHw0yziJ2upIYWAAUhAA7kAQgo2bf1z++brxlsGl9wgTWuCYB4C/wEkBbNlXl5SPzxlT8KVumD4GAGqA8wLYrG0p7OnDd26c/2d/TB8DADVAEIJtOpXDpu2RJIWIqK8ftzpSxPwxAFAj09xOfObMmePHjzPGhg4dGhERUb1AfHx8VlaWYdnd3X3MmDEm2S5AjdKK2cQ4uUhPRNTanduCG+cBoHYmOD1888037733XlRUlCAI8+fP//jjj+fOnVulzKJFi0pLS8PCwogoJCQEQQjNp0BPY3fIqcWMiDwdaVOMEOisdp0AwIKZIAgnTpw4Z84cQRCIaMiQIQsWLKgehET00ksvTZkypembA7gLvUJTdkkJeYyIHHj6OUrs6o0uUQC4GxMEoZ+fX8Wyu7s7x9V83jl9+rQkSREREd27d2/6RgFq9ORBOS6NERFHtGKYMCoEKQgAdTDllROdTrdw4cJnnnmm+lu+vr6XL19OSkp69tlnJ06c+M0339T2IYqi/OMf/6h42aNHj+jo6Lts0cHBoYnVhsbR6/WKoqhdi79Ycp779oox+d7uwaaE63U6dWvUjHQ6naOjo9q1sFPY+SrS6XSyLMuyXM/yDg4OPF/HsND6BmHLli2rr3zvvfdmzZplWJYkaebMmS1atHjxxRerl1y7dq1hIT09vWvXrtOnT4+JialtW7dv365YliSpnjUEO7cphd4+Y0zB2W3Z692ZuvUBAGvBMVav80VaWlr1lV5eXq6urkQky/KcOXNycnI2btyo0dQxpX9MTMzYsWNfeOGF6m9JkuTs7KzX6+tTJSIqKipyc3OrZ2EwLUOLsM5/bvM4k8uGbJaKJSKioUFcXKyosfWbBgsLC93d3dWuhZ3CzleRVqsVRVEUTdmdWd/PCg0Nre0txtjTTz+dkZGxefPmyqfFvLw8nU4XHBzMGKu4cFhQUHDu3Llnn322KZUGqCyjhO7dKRtSsLU793OU7acgAJiQCUL1yy+//PLLL8eNG/fwww8b1qxcudLBweHjjz9OSEjYvHlzenr6+PHjhw8f7uDgsHHjxh49eowbN67p2wUgolKJ7tslpRQzIvJwoM0xgr+T2nUCAKtS367Ru7h8+XJCQkLlNffffz/P84mJifn5+YMGDZIkae/evYmJiUQUERERHR1d28hSdI1aEUvoGmVED+2VV11TiEjkadtoMSrUXoaJondORdj5KmqOrlETBKEJIQitiCUE4Vsn5fdOG0eufj5YeKqzHU0ZiHOxirDzVdQcQWhHJw6wMWuTlH+Up+ALXXm7SkEAMCGcO8AqHc5kD8fLht6M0WHcR3iyBAA0FoIQrE9yEZu8S9LKRESdvbg1kaKIAxkAGgvnD7AyhXoav0POLCUiCnCmraMFT0zxAQBNgCAEa6IwenCvdP4WIyIngTZEia3c7WWYKAA0EwQhWJPXfpe33DCOc/5qqDAoECkIAE2FIASr8cNV5aME4zDR13vwM9vh6AUAE8CpBKzD4Uz2+EHjfPOx4dzivhgmCgCmgSAEK3CjiE3eJelkIqIIL27VSFFAnygAmAiCECxdqUT37zIOE/XV0KYYDBMFAFNCEIJFY0Rz9ssnchgROfC0Pkps64HGIACYEoIQLNpbJ+S1ScYBMp8NEkYEIwUBwMQQhGC51l1X3j9jTMGXu/GPdcLhCgCmhzMLWKhTOWxOpdlE/w+ziQJA80AQgiXKKKGJcXKJRETUyYtbHYlhogDQXBCEYHFKJZoUJ6UWMyLy0dDmGMELw0QBoNkgCMGyMKJHD8jHs43DRNeNEtthmCgANCcEIViWRaeUVdeMA2Q+HShEhiAFAaB5IQjBgvzyp7LwlHEetee68E/iofMA0PxwogFLcTqXzd5nHCYaHcp9MgDDRAHAHBCEYBEySujenXKxRETUwRMPnQcA88HJBtRXItHE8mGi3hraHCN4a9SuEwDYDQQhqExh9HC8/Hs2IyKRpzWRYgdPDJABAPNBEILK3jwhr79uHCb673uE6FCkIACYFYIQ1PT9VeWDs8YUfLEr/3QEDkgAMDecd0A1B26yJyo9dP4jDBMFADUgCEEdSYXs/vKHznfxxkPnAUA1CEJQQZ6Oxm6Xs7VEREHOtHU0HjoPAKpBEIK56RWatlu6nM+IyFmkjdFiCzc0BgFANQhCMLfnDsu70xkRcUTfDhUGBCAFAUBNCEIwq38mKF9cMg4Tfb+fMKMtjkAAUBlOQ2A+v6WwN343DhN9uD3/Wg8cfgCgPpyJwExO57LpuyXDpNpDg7gvhuBmCQCwCAhCMIf0ElYxp3Ybd+7nKFGDHAQAy4AghGZXKtGkONkwp7aPhraNEfyd1K4TAEA5BCE0L4XRg3uNc2o78LR2FObUBgDLgiCE5jX/uLwx2ThM9D+DhFEhSEEAsCwIQmhG315RPj5nTMFXe/BPdMLxBgAWBycmaC77b7KnyufUntyKf78vhscAgCVCEEKzuHSbTYqTyhQiot5+3A8jBB59ogBgkRCEYHq5Oro3Tr6lIyIKceF+jRZcRbXrBABQCwQhmJheoam7pKv5jIhcRNoYLYS5ojEIAJYLQQimxIge3S/vzWBExHP0v5FCP3+kIABYNAQhmNI7J+Uf/zAOE/2ovzCxJQ4wALB0OE+ByXx1SVl82piCj3TgX+qGowsArABOVWAaW26wpw8bb5YYG445tQHAaiAIwQSOZ7MZeyRJISLq68etiRRFHFkAYCVwuoKm+qOAjd8hGZ4s0daD2zJadHNQu04AAPWGIIQmSS9h0dvkbC0Rkb8TbRstBDqrXScAgIbAfc7QeAV6Ghcn/1nIiMjdgbaPEdvjyRIAYG3QIoRGKlPogXjuTK7x+UrrRom9/ZCCAGB9EITQGIzoiUO0O4MjIo7oyyHC6DCkIABYJQQhNMbzR+SVScblD/oJczrgQAIAa4XzFzTYWyfl/yQab5z/exf+1R44igDAiuEUBg2z5JzyXvn0MdNa09J7cOM8AFg3BCE0wHdXlPnHjNPHRIfQVwMVPGUQAKwdghDq65c/lbkHZEZERIMCubUjSYPWIABYPwQh1MvONPbgXtkQg/38ue1jRDxrFwBsA4IQ6rY3g03cKelkIqKu3ty2MaI7JlEDAFuBIIQ6HM9mE3dKWpmIqK0HtyNW8NWoXScAANNBEMLdnL/Fxm6XCvVERKGuXFysEOKC4TEAYFMQhFCrC7dZ1FYpV0dE5O9EcbFCa3ekIADYGgQh1Oz8LTbyNymzlIjI05G2jxE7eyEFAcAGYeQf1OBcHhu1VTI8XMnDgbaNwYTaAGCz0CKEqs7mscjyFPR0pB2x4sAApCAA2Cy0COEvzuSy6G1STkUKjhEHIAUBwKYhCOGO07ksunx0jJcj7YgV+/sjBQHAxiEIwahKCu6MFfshBQHADiAIgYjoVA6L3ibl6YiIvDW0M1bsi9ExAGAfMFgG6CRSEADsmAlahDdu3Ni+fXvFy9jY2PDw8CplZFn+4YcfLl26FBERMXPmTEHAYwssxd4MNmmnVKAnIvJ3ol1jxe4+SEEAsCMmaBEmJCS8/fbbJ8vl5+dXL/P4448vX748PDz8888/f+qpp5q+UTCJddeV2O3GFAxwpj3jkIIAYHdMc42wbdu2X3zxRW3vpqSkrFy5Mjk5OTAwcMqUKa1bt37nnXdCQ0NNsmlotGUXlOePyAojIgpz5XbEChGYOwYA7I9prhHm5OR88sknK1asuHnzZvV3Dx482LVr18DAQCIKCgqKiIg4ePCgSbYLjcOI3j0lP3fYmIKdvbhDE5CCAGCnTNAidHFx6d69e1ZW1oEDB1566aWtW7cOHDiwcoGMjAx/f/+Kl4GBgenp6bV9mqIoc+fOrXg5ZMiQGTNm1FZYq9WKIga+NozM6Pnj/DdXjbHX15dtjFR8RUmrbdjn6PV6RVEYY6avItSDVqt1cMBjIdWBna8iw2m//md+BweHOkel1Ouz4uPjo6Kiqq8/e/ZsREREZGRkZGSkYc0bb7zxxhtv7N279y/bEEVZlite6vX6uxxDHMf17du34mX79u3vUtjBwQGHY4PoZJq9n/3ypzG9xodz/xvBuzT2bwlFUbD/1YKDX0XY+SqSZblBQcjzdXd81uuzhg8fXlZWVn09x1XtTBsyZMjKlSurrAwJCancBExPT7/LBUKO45588sn61IqIBEHAANT6u6Wje+OkgzeNKTi7Pf/1UMGhsb3jiqJwHIf9rxYc/CrCzleRUM6En1nfsyBXE8NbOp2uotjWrVsjIiIMy2fPnk1NTSWiUaNGXb9+/eLFi0R04cKF5OTkkSNHmvB3gPpIL2EjfruTgn/vwn83vPEpCABgM0xwge1vf/tbZmZmixYtLl68mJqaWnFP4dNPP33ffffNmzfP29t7wYIFo0ePjo2N3bZt21tvveXl5dX07UL9JeSxcTvk1GJGRDxH/x4oPBuBDAQAICLimj7Y4datW0ePHs3Ozg4ODh48eLCLi4th/fnz5318fEJCQgwvT506lZiY2LVr1169etX2UZIkOTs76/X6em66qKjIzc2tifW3eZuSlYf2yUV6IiKNQD+OEKa2NkEKGgbLaDSapn8UNEJhYaG7u7vatbBT2PkqauhgmfowQRCaEILQ5D5KUF773XibhKcjbYgWRwab5jYJBKG6cC5WEXa+ipojCHHvgc3SyfTUIXnFFcXwsq0HtykGNwsCAFSFILRNacXs/t3ysSxjc39YEPdzlOjnpG6lAAAsEYLQBh28yabulm6WGl8+0oH/7xDBEYNjAABqgiC0Nf9JVF4+JusVIiIHnj4eIDzXBRkIAFArBKHtKNDT3P3yuuvGi4KBzrR2lDgsCBcFAQDuBkFoI87msam75av5xouC/f25n6OEMFekIABAHRCEtuC/F5WXj8mlkvHlMxH8xwMEDWaAAgCoBwShdcvR0qMH5E3Jxu5Qdwf6aqgwvQ0uCgIA1BeC0IrFpbE58XJ6ibE7tIcPt3aU0MET3aEAAA2AILRKBXp664S87IJimDKGI3qsE7/0HqHRD1QCALBbOHFan99S2JMHjTNoE1GgM303XBwThoYgAEBjIAitSY6WXjgqr/xDqVhzXyt++WAhwFnFSgEAWDcEodXYfEN58qBScUUwwJk+6i/Mbo9xMQAATYIgtAIZJfT0IXlj8p2G4Mx2/L8GCr548AMAQJMhCC1amUL/vaC8c0rOLzOuCXahzwcLk1qiIQgAYBoIQgvFiNYmKa/9rvxZaOwL5Yge78T/c4Dg4aBu1QAAbAqC0BIdzWIvHZWPZN15ZnKEF/fZYGGEiZ6pCwAAFRCEliW1mL3xu/LTH0pFBvpo6O1ewjMRvIjeUACAZoAgtBTZWlqSIH+aqGhl4xqNQH/vwi/oKXg6qlozAACbhiBUX2YpfZQgL7+oFJfPms0RTWnN/19/vrU7+kIBAJoXglBNacVsyTnly0tKiXRnZX9/7pN7hMGBiEAAAHNAEKrjZA5bel5Zm6To79wcSL18ubd68ZNa8chAAACzQRCalaTQ5hvK0vPKgZus8vq+ftzbvfnxLRCBAADmhiA0k6RC9u1lZcUVVjFHmsGIYO6V7sLYcCQgAIA6EITNq0SiX5OVby4re9JZ5QB05Gl6G/7FbnwvX0QgAICaEITNQifTjlRldRLblHxnLKhBkDM92pF/JkIIdlGpcgAAUAmC0JRu6SguTfkthW1KVm6X/eUtgaMxYdzcTvz4cNwaDwBgQRCEJnA2j21LYVtTlCNZTFKqvtvFm3ugLf9wey7MFb2gAAAWB0HYSJfz2ZFMduAm25HG0opZ9QJtPbjpbbgZbfhuPsg/AADLhSCsL61MJ7LZoUx2OIsdyVSytTWU4Tnq58+NDefHhXN9/JB/AABWAEFYs1s6SipkibfYhdssqYASb7ErBTV0exr4aGhUCB8Vyo1vwYW4IP8AAKwJgpCKJbqSz67ks8v5dPk2u5zPruSzQn0dP+XnRAMD+MGB3NAgbkAAJyD+AACsk10EYbFE2aXsZinlaClHyzJLKUvLskspvYRdyaeUmq7wVcdz1NGTGxTIDQ7kBgVyHT0RfQAAtsC6g7BUotwidruMbuvoVhm7paM8Hd0sYdlaytGybC3dLKUcLSuR6v6oKrwcqaMX19GT6+jJdfCkDp5cB0/OSWiG3wEAAFRlxUH4v2vK3+I1ZUrDU+6vHHhq48519OI6elIHT2P4BTibpI4AAGDprDgIV19jZbWMXqnCSSB/Jy7AmQKcyd+J83OiQGcuwIn8nbn2HtTanXPAHe4AAPbKioPw+a782VyZON7LkbwcyVvDeTmSj4b8nblAZ/J3Ij8nLsCJAp05Nwe16woAAJbKioNwVAiXOKHMzc1N7YoAAIAVQ58gAADYNQQhAADYNesOwh9//DE9PV3tWtip48eP7969W+1a2KmCgoLly5erXQv7tWzZstLSUrVrYae2bdt25swZ036mdQfhmjVrLl26pHYt7NTBgwd37dqldi3sVGpq6rfffqt2LezX8uXLMzMz1a6FndqxY8fRo0dN+5nWHYQAAABNhCAEAAC7hiAEAAC7xjFWrymnzUOWZY1G07Jly3qWz8jI8Pb2dnJyatZaQY3y8/NlWfbx8VG7IvZIr9dnZmaGhYWpXRE7lZKSEhISIgiYfVgFubm5jo6O7u7u9Sz/4IMPLl68+O5lLOuGekEQkpKSJKm+04fqdDqNRtOsVYLayLLMGBNFyzqE7AcOfhVh56tIkiSO4+r/V0hwcHCdZSyrRQgAAGBmuEYIAAB2DUEIAAB2DUEIAAB2DUEIAAB2zWqG/CUnJ69YsaKkpGTatGl9+/atXkCr1X799ddJSUl9+/Z94IEHOI4zfyVtkk6n27Vr17FjxyRJGjZs2JgxY6qX2b59+40bNwzLbm5uDz74oHnraMtWrVpVWFhoWA4PD4+Nja1e5tSpU2vWrHFycnr44YfbtGlj3grasg0bNmRnZ1e8DAoKuvfeeysXyMvLW79+fcXLIUOGREREmK9+NkdRlMuXL58+fbqoqOjxxx+v/NauXbu2bt0aEBAwd+5cPz+/6j+blJT0/fffa7XaGTNm9OrVq0HbtY4WYUZGRt++fQsLCwMCAqKiog4cOFC9zP33379p06b27dt/+OGHb7zxhvkraavWrFmzePFiQRA8PT3nzp27YMGC6mWWLVu2efPmpKSkpKSkikQEk1iwYMHhw4cN+/bmzZvVCxw9enTkyJE+Pj5arbZ///6pqanmr6StSktLSyr34Ycfbt26tUqB1NTUefPmVZQpKChQpZ424/Dhw5GRkZ9//vkTTzxRef3KlStnzZrVqlWrS5cuDRo0qPqM5ykpKf369dNqtT4+PiNGjDh27FjDNsyswbvvvjt58mTD8ocffjh+/PgqBU6dOuXh4VFcXMwYu3DhgpubW35+vrlraaNKS0srlrdv3+7l5VW9zLhx43744QczVsqOtG7d+sSJE3cpMHny5IULFxqWH3zwwTfeeMMs9bIvhYWFbm5ux48fr7L+7NmzYWFhqlTJJhnuTr548WKVbOrSpcvq1asZY4qi9OnT5/vvv6/yg6+//voDDzxgWF64cOH999/foO1aR4swPj4+JibGsBwdHR0fH1+9wODBg11cXIioc+fO3t7eJ0+eNHctbVTliXu0Wq2bm1uNxfbu3btkyZJff/1VlmVzVc1erFu3bunSpXv37q3x3fj4+OjoaMNyjd8OaLpVq1a1atWqX79+1d8qLS1dunTp8uXLr169av6K2RieryGScnJyEhMTDQc5x3FRUVHVD/L9+/ffPSPq2G6jamtuGRkZ/v7+huWAgIDCwsKioqLKBW7evFlRwFAGzyk0uaKiotdee+3VV1+t/lb79u3d3NxycnJef/31yMhIvV5v/urZqr59+8qynJKS8tBDDz322GNV3tXpdHl5eZW/HRkZGWavo+375ptvHn300errHR0dhwwZkpube+TIkd69e//yyy/mr5vNy8jIEEXR29vb8DIwMLD66b1KRuTm5paVldV/E9YxWMbBwaFi3jXD/DpVZvYSRbFyQ0Sv1zs6Opq1iraurKxs2rRp/fr1e/rpp6u/u3TpUsPCW2+9FRER8fPPP8+YMcO8FbRZa9euNSy88MILHTp0+Pvf/96tW7eKd0VR5Hm+8rcDR77JXbp06fTp05s3b67+VqdOnTZu3GhYHjx48Lx58yZPnmze2tk+BwcHRVEYY4YhkHq9vvr8dqIoVv4W8DzfoJlgraNFGBoaWvEnQFpamq+vb5WJtkNDQ9PS0gzLjLH09PSQkBBz19J2lZWVTZ482dPTc8WKFTX2XVRwdXXt2bPn9evXzVY3+9GiRYuwsLAq+1YQhICAgIqDPy0tDUe+yX311VcTJ06s3OdUo0GDBt24cQOXBkwuJCSEMVbR1ZGWllZ9+tAqGREUFGSDQThhwoT169crikJE69atmzBhgmH90aNHDePoxo4de+zYMcPpYN++faIo9u/fX8UK2xK9Xj9t2jSNRvPjjz9WPraSk5PPnDlDRIqiVPRCZGdnHz16tEuXLurU1ebodDrDYU9ECQkJN27cMIzOz8rKOnz4sGH9vffeaxjBzxhbv359xbcDTKKsrOynn3565JFHKq+Mj4+/desWEWm12oqVmzdv7tSpEx5JYXIeHh7Dhg1bt24dEZWWlm7ZssVwE0tRUdGePXsMX5AJEyasW7eOMUZ/zYj6atTQHnMrLCzs1avX8OHDp02bFhQUdOXKFcP6du3aVQwfevnll9u0afO3v/0tMDBwxYoVqtXV5nzxxRdE1L179z7lCgoKGGOLFy8eNWoUYyw3N9ff33/SpEkzZszw8/ObOXOmoR8Dmm7//v2tW7eeOnXqfffd5+7uvmjRIsP6VatWtWjRwrB87dq14ODgKVOmREZGduvWDeOlTWv9+vVhYWGSJFVe6ezsvGvXLsbYK6+8MmDAgJkzZw4dOtTPzy8+Pl6latqI/Pz8Pn36GP6S7tOnT3R0tGH9gQMHfH19Z82a1bt379jYWMPg0tOnTxuikTF2+/btrl27RkZGTpkyJSQkJCkpqUHbtZqnT+h0ut27dxcXF0dFRVVcNU1ISAgNDfX19TW8PH78+LVr13r37t2xY0f1amprsrKyUlJSKq/p2bOnIAjp6emFhYWGXX316tXz589LkhQREYHmoAkpinLu3LnLly87Ojr26tWr4lGdeXl5KSkpPXr0MLy8ffv27t27nZycRo0ahcdzmlZqampZWVmVaQpOnjzZvn17Dw+PkpKS48ePZ2Rk+Pn5DRgwwMPDQ6162gZZlg39TAYODg7du3c3LGdkZOzfv9/Pz2/kyJGGCzQlJSXnz5/v16+f4dqhVqvdvXu3VquNiory9PRs0HatJggBAACag3VcIwQAAGgmCEIAALBrCEIAALBrCEIAALBrCEIAALBrCEIAALBrCEIAALBrCEIAALBrCEIAALBrCEIAALBrCEIAALBr/w+77sbW81iNOQAAAABJRU5ErkJggg==",
"image/svg+xml": [
"\n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n",
" \n",
" \n",
" \n",
"\n",
"\n",
" \n"
],
"text/html": [
"Convergence Theorem (Relaxation) \n",
"\n",
"Suppose that $f: [a,b]\\to\\mathbb R$ is continuously differentiable with $f(\\xi) = 0$, $f'(\\xi) \\not= 0$. Then, there exists $\\delta, \\lambda> 0$ such that $x_{n+1} = x_n - \\lambda f(x_n)$ converges at least linearly to $\\xi$ for all $x_1 \\in I_\\delta := [\\xi-\\delta,\\xi+\\delta]$. The asymptotic error constant (for $\\alpha =1$) is given by $\\left| 1 - \\lambda f'(\\xi) \\right|$.\n",
"\n",
"
\n",
"\n",
" Example (Kepler's equation). \n",
"\n",
"\n",
"We have,\n",
"\n",
"\\begin{align*}\n",
" &f(\\psi) = \\psi - \\epsilon \\sin \\psi - 2\\pi \\qquad &f(2\\pi) = 0 \\nonumber\\\\\n",
" &f'(\\psi) = 1 - \\epsilon \\cos \\psi \\qquad &f'(2\\pi) = 1 - \\epsilon. \\nonumber\\\\\n",
"\\end{align*}\n",
"\n",
"Therefore, for $\\epsilon \\not= 1$, we have $f(2\\pi) = 0 \\not= f'(2\\pi)$ and so we can apply the above theorem. Relaxation converges linearly with asymptotic error constant $\\mu := \\left| 1 - \\lambda f'(2\\pi) \\right| = \\left| 1 - \\lambda (1-\\epsilon) \\right|$. \n",
"\n",
"In the above example, we chose $\\epsilon = 0.9$ and so the asymptotic error constant is $\\left| 1 - 0.1\\lambda \\right|$. This is in good agreement with what we see numerically.\n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "aa1ab311",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"┌───────────┬────────────────┬───────────┬─────────┬────────────┐\n",
"│\u001b[1m iteration \u001b[0m│\u001b[1m absolute error \u001b[0m│\u001b[1m log error \u001b[0m│\u001b[1m alpha \u001b[0m│\u001b[1m mu (α = 1) \u001b[0m│\n",
"├───────────┼────────────────┼───────────┼─────────┼────────────┤\n",
"│ 1.0 │ 0.283185 │ -0.547929 │ NaN │ NaN │\n",
"│ 2.0 │ 0.251474 │ -0.599507 │ 1.09413 │ 0.888019 │\n",
"│ 3.0 │ 0.223949 │ -0.649852 │ 1.08398 │ 0.890544 │\n",
"│ 4.0 │ 0.199873 │ -0.699245 │ 1.07601 │ 0.892496 │\n",
"│ 5.0 │ 0.178691 │ -0.747898 │ 1.06958 │ 0.89402 │\n",
"│ 6.0 │ 0.159967 │ -0.795969 │ 1.06427 │ 0.895218 │\n",
"│ 7.0 │ 0.143357 │ -0.843581 │ 1.05982 │ 0.896166 │\n",
"│ 8.0 │ 0.12858 │ -0.890827 │ 1.05601 │ 0.89692 │\n",
"│ 9.0 │ 0.115403 │ -0.937782 │ 1.05271 │ 0.897522 │\n",
"│ 10.0 │ 0.103633 │ -0.984504 │ 1.04982 │ 0.898004 │\n",
"│ 11.0 │ 0.0931025 │ -1.03104 │ 1.04727 │ 0.89839 │\n",
"│ 12.0 │ 0.0836712 │ -1.07742 │ 1.04499 │ 0.8987 │\n",
"│ 13.0 │ 0.0752163 │ -1.12369 │ 1.04294 │ 0.89895 │\n",
"│ 14.0 │ 0.0676308 │ -1.16986 │ 1.04109 │ 0.899152 │\n",
"│ 15.0 │ 0.0608214 │ -1.21594 │ 1.0394 │ 0.899314 │\n",
"│ 16.0 │ 0.0547055 │ -1.26197 │ 1.03785 │ 0.899445 │\n",
"│ ⋮ │ ⋮ │ ⋮ │ ⋮ │ ⋮ │\n",
"└───────────┴────────────────┴───────────┴─────────┴────────────┘\n",
"\u001b[36m 14 rows omitted\u001b[0m\n"
]
}
],
"source": [
"orderOfConvergence( x, 2π; α=1 )"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "6672cde9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f'(ξ) = 0.09999999999999998\n",
"Theoretical asymptotic error constant μ = 1 - λf'(ξ) = 0.9\n"
]
}
],
"source": [
"println( \"f'(ξ) = \", f_prime(ξ))\n",
"println( \"Theoretical asymptotic error constant μ = 1 - λf'(ξ) = \", 1 - λ*f_prime(ξ) )"
]
},
{
"cell_type": "markdown",
"id": "c6b67325",
"metadata": {},
"source": [
" Proof. \n",
" \n",
"Relaxation is a simple iteration with $g(x) = x - \\lambda f(x)$. We wish to apply the contraction mapping theorem and thus it is sufficient to check $|g'(\\xi)| \\leq L < 1$.\n",
"\n",
" Suppose that $f'(\\xi) > 0$ (the proof is very similar in the other case). Since $f'$ is continuous, there exists $\\delta> 0$ such that \n",
"\n",
" \\begin{align*}\n",
" m := \\frac12 f'(\\xi) \\leq f'(x) \\leq M := \\max_{y\\in [a,b]} f'(y).\n",
" \\end{align*}\n",
"\n",
"As a result, we have $1 - \\lambda M \\leq 1 - \\lambda f'(x) \\leq 1 - \\lambda m$. Choosing $L = 1 - \\lambda m = \\lambda M - 1$ gives $\\lambda = \\frac{2}{M+m}$ and thus \n",
"\n",
"\\begin{align*}\n",
" L = 1 - \\frac{2m}{m + M} = \\frac{M-m}{M+m} < 1.\n",
"\\end{align*}\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "ee25d30e",
"metadata": {},
"source": [
"## Newton\n",
"\n",
"* $x_{n+1} = x_n - \\frac{f(x_n)}{f'(x_{n})}$"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "5959cd06",
"metadata": {
"jupyter": {
"source_hidden": true
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to c:\\Users\\math5\\Math 5485\\Lectures 1-12\\Pictures\\Newton.gif\n",
"└ @ Plots C:\\Users\\math5\\.julia\\packages\\Plots\\xKhUG\\src\\animation.jl:156\n"
]
},
{
"data": {
"text/html": [
"Theorem. \n",
"\n",
"Suppose $f: [a,b] \\to \\mathbb R$ is twice continuously differentiable with $f(\\xi) = 0$ and $f'(\\xi) \\not= 0$ for some $\\xi \\in [a,b]$. Further suppose that \n",
"\n",
"\\begin{align}\n",
" \\left|\\frac{f''(x)}{f'(y)}\\right| \\leq A\n",
"\\end{align}\n",
"\n",
"for all $x,y \\in [a,b]$. Then, the Newton iteration $x_{n+1} = x_n - \\frac{f(x_n)}{f'(x_n)}$ converges at least quadratically to $\\xi$ for all $x_1 \\in [a,b]$ such that $|x_1 - \\xi| \\leq A^{-1}$. The asymptotic error constant $\\mu$ is $\\frac12 \\left|\\frac{f''(\\xi)}{f'(\\xi)}\\right|$.\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "226ceecb",
"metadata": {},
"source": [
"Example. \n",
"\n",
"\n",
"Notice that you have already seen an example of Newton's method! \n",
"\n",
"\\begin{align}\n",
" x_{n+1} = \\frac12 \\left( x_n + \\frac{a}{x_n} \\right)\n",
"\\end{align}\n",
"\n",
"which is the Newton iteration with $f(x) = x^2 - a$. Notice that $\\xi = \\sqrt{a}$ and $f'(\\xi) = 2\\sqrt{a}$, $f''(\\xi) = 2$. You showed that the order of convergence is quadratic in this case. Compare with the theorem\n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "7ee07def",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"┌───────────┬────────────────┬───────────┬─────────┬────────────┐\n",
"│\u001b[1m iteration \u001b[0m│\u001b[1m absolute error \u001b[0m│\u001b[1m log error \u001b[0m│\u001b[1m alpha \u001b[0m│\u001b[1m mu (α = 2) \u001b[0m│\n",
"├───────────┼────────────────┼───────────┼─────────┼────────────┤\n",
"│ 1.0 │ 0.414214 │ -0.382776 │ NaN │ NaN │\n",
"│ 2.0 │ 0.0857864 │ -1.06658 │ 2.78644 │ 0.5 │\n",
"│ 3.0 │ 0.0024531 │ -2.61028 │ 2.44734 │ 0.333333 │\n",
"│ 4.0 │ 2.1239e-6 │ -5.67287 │ 2.17328 │ 0.352941 │\n",
"│ 5.0 │ 1.59472e-12 │ -11.7973 │ 2.0796 │ 0.353522 │\n",
"└───────────┴────────────────┴───────────┴─────────┴────────────┘\n",
"theoretical value of μ (for α=2) = 0.35355339059327373\n"
]
}
],
"source": [
"g(x) = x^2 - 2\n",
"dg(x) = 2*x \n",
"orderOfConvergence( Newton(g,dg,1.), sqrt(2); α=2 )\n",
"\n",
"println( \"theoretical value of μ (for α=2) = \", (1/2) * ( 2/ dg(sqrt(2)) ) )"
]
},
{
"cell_type": "markdown",
"id": "071a038f",
"metadata": {},
"source": [
"Proof. \n",
"\n",
"First, notice that $x_{n+1} = x_n - \\frac{f(x_n)}{f'(x_{n})}$ and so \n",
"\n",
"\\begin{align}\n",
" -f(x_n) &= f'(x_{n}) (x_{n+1} - x_{n}) \\nonumber\\\\\n",
" &= f'(x_{n}) (\\xi - x_n) + f'(x_{n}) (x_{n+1} - \\xi). \n",
"\\end{align}\n",
"\n",
"At the same time, there exists $\\eta_n$ between $x_n$ and $\\xi$ such that \n",
"\n",
"\\begin{align}\n",
" 0 = f(\\xi) \n",
" &= f(x_n) + f'(x_n) ( \\xi - x_{n} ) + \\frac{1}{2} f''(\\eta_n) (x_n - \\xi)^2 \\nonumber\\\\\n",
" &= f(x_n) + \\big[ -f(x_n) - f'(x_{n}) (x_{n+1} - \\xi) \\big] + \\frac{1}{2} f''(\\eta_n) (x_n - \\xi)^2.\n",
"\\end{align}\n",
"\n",
"Therefore, assuming $x_n \\in [a,b]$, we obtain \n",
"\n",
"\\begin{align}\n",
" \\left| x_{n+1} - \\xi \\right| &= \\frac{1}{2} \\left| \\frac{f''(\\eta_n)}{f'(x_n)} \\right| (x_n - \\xi)^2 \\nonumber\\\\\n",
" %\n",
" &\\leq \\frac12 A |x_n - \\xi|^2.\n",
"\\end{align}\n",
"\n",
"If $|x_1 - \\xi| \\leq A^{-1}$ then, we have $|x_2 - \\xi| \\leq A |x_1 - \\xi|^2 \\leq A A^{-1} |x_1 - \\xi| \\leq A^{-1}$. Proceeding iteratively, we obtain $|x_{n} - \\xi| \\leq A^{-1}$. As a result, we have $\\left| x_{n+1} - \\xi \\right| \\leq \\frac12 |x_n - \\xi|$ and so $(x_n) \\to \\xi$. Finally, we have \n",
"\n",
"\\begin{align}\n",
" \\frac{\\left| x_{n+1} - \\xi \\right|}{|x_n - \\xi|^2} = \\frac{1}{2} \\left| \\frac{f''(\\eta_n)}{f'(x_n)} \\right| \\to \\frac12 \\left| \\frac{f''(\\xi)}{f'(\\xi)} \\right|\n",
"\\end{align}\n",
"\n",
"as $n\\to\\infty$.\n",
"
\n",
"\n",
" Example (Kepler's equation). \n",
"\n",
"Recall that\n",
"\n",
"\\begin{align*}\n",
" &f(\\psi) = \\psi - \\epsilon \\sin \\psi - 2\\pi \\qquad &f(2\\pi) = 0 \\nonumber\\\\\n",
" &f'(\\psi) = 1 - \\epsilon \\cos \\psi \\qquad &f'(2\\pi) = 1 - \\epsilon \\nonumber\\\\\n",
" &f''(\\psi) = \\epsilon \\sin \\psi \\qquad &f''(2\\pi) = 0 \\nonumber\\\\\n",
" &f'''(\\psi) = \\epsilon \\cos \\psi \\qquad &f'''(2\\pi) = \\epsilon.\n",
"\\end{align*}\n",
"\n",
"Therefore, for $\\epsilon \\not= 1$, we can apply the above theorem. Newton therefore converges at least quadratically with asymptotic error constant $\\mu = \\frac12 \\left| \\frac{f''(2\\pi )}{f'(2\\pi)} \\right| = 0$. As a result, we would expect Newton iteration to converge super-quadratically in this case:\n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "cea26522",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"┌───────────┬────────────────┬───────────┬─────────┬────────────┐\n",
"│\u001b[1m iteration \u001b[0m│\u001b[1m absolute error \u001b[0m│\u001b[1m log error \u001b[0m│\u001b[1m alpha \u001b[0m│\u001b[1m mu (α = 2) \u001b[0m│\n",
"├───────────┼────────────────┼───────────┼─────────┼────────────┤\n",
"│ 1.0 │ 0.783185 │ -0.106135 │ NaN │ NaN │\n",
"│ 2.0 │ 0.374019 │ -0.427107 │ 4.02417 │ 0.609767 │\n",
"│ 3.0 │ 0.0954133 │ -1.02039 │ 2.38908 │ 0.68206 │\n",
"│ 4.0 │ 0.00250109 │ -2.60187 │ 2.54988 │ 0.274733 │\n",
"│ 5.0 │ 4.69349e-8 │ -7.3285 │ 2.81663 │ 0.00750305 │\n",
"└───────────┴────────────────┴───────────┴─────────┴────────────┘\n"
]
}
],
"source": [
"orderOfConvergence( y, 2π ; α = 2)\n",
"# plot( abs.( y .- ξ ) , yaxis=:log, xaxis=:log, m=:o, lw=3, title=\"absolute error\", legend=false )"
]
},
{
"cell_type": "markdown",
"id": "4ebf1a22",
"metadata": {},
"source": [
" \n",
"\n",
"In fact, in this case the convergence is cubic with asymptotic error constant $\\mu = \\frac13 \\left| \\frac{f'''(\\xi)}{f'(\\xi)} \\right| = \\frac13 \\frac{\\epsilon}{1-\\epsilon}$ (**Exercise**: why?). \n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "8d8526c5",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"theoretical value of μ (for α=3) = 3.0000000000000004\n"
]
}
],
"source": [
"println( \"theoretical value of μ (for α=3) = \", (1/3) * ( ϵ/ (1-ϵ) ) ) "
]
},
{
"cell_type": "markdown",
"id": "7d0ded19",
"metadata": {},
"source": [
"\n",
"\n",
"In your assignment, you will experiment with the case $\\epsilon = 1$ and so $f(2\\pi) = f'(2\\pi) = f''(2\\pi) = 0$.\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "27134deb",
"metadata": {},
"source": [
"## Secant Method\n",
"\n",
"* Consider $x_{n+1} = g(x_n, x_{n-1})$ \n",
"* Consider Newton but with an approximate derivative: \n",
"\n",
"\\begin{align}\n",
" f'(x_n) &\\approx \\frac{f(x_{n}) - f(x_{n-1})}{x_{n} - x_{n-1}}\n",
" \\nonumber\\\\ &= f[x_n, x_{n-1}] \n",
"\\end{align}\n",
"\n",
"* That is, $x_{n+1} = x_n - \\frac{f(x_n)}{f[x_n, x_{n-1}]}$"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "b79de49b",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"┌ Info: Saved animation to c:\\Users\\math5\\Math 5485\\Lectures 1-12\\Pictures\\Secant.gif\n",
"└ @ Plots C:\\Users\\math5\\.julia\\packages\\Plots\\xKhUG\\src\\animation.jl:156\n"
]
},
{
"data": {
"text/html": [
"Theorem. \n",
"\n",
"Suppose $f: [a,b] \\to \\mathbb R$ is continuously differentiable with $f(\\xi) = 0$ and $f'(\\xi) \\not= 0$ for some $\\xi \\in [a,b]$. Then, the Secant iteration $x_{n+1} = x_n - \\frac{f(x_n)}{f[x_n,x_{n-1}]}$ converges at least linearly to $\\xi$ for all $x_1,x_2$ sufficiently close to $\\xi$. \n",
"
\n",
"\n",
"Remark. \n",
"\n",
"In fact, order of convergence $\\alpha = \\frac{1 + \\sqrt{5}}{2}$ in this case, with $\\mu = \\left| \\frac{f''(\\xi)}{2 f'(\\xi)} \\right|^{\\alpha/(1 + \\alpha)}$. \n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "c28af500",
"metadata": {},
"source": [
"Example. \n",
"\n",
"Let's try and approximate $\\sqrt{2}$ again:\n",
"\n",
"
"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "f004616f",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"┌───────────┬────────────────┬───────────┬─────────┬────────────────────────────\n",
"│\u001b[1m iteration \u001b[0m│\u001b[1m absolute error \u001b[0m│\u001b[1m log error \u001b[0m│\u001b[1m alpha \u001b[0m│\u001b[1m mu (α = 1.618033988749895\u001b[0m ⋯\n",
"├───────────┼────────────────┼───────────┼─────────┼────────────────────────────\n",
"│ 1.0 │ 0.414214 │ -0.382776 │ NaN │ Na ⋯\n",
"│ 2.0 │ 0.585786 │ -0.232261 │ 0.60678 │ 2.4382 ⋯\n",
"│ 3.0 │ 0.0808802 │ -1.09216 │ 4.70229 │ 0.19215 ⋯\n",
"│ 4.0 │ 0.0142136 │ -1.8473 │ 1.69142 │ 0.83147 ⋯\n",
"│ 5.0 │ 0.000420584 │ -3.37615 │ 1.82761 │ 0.41005 ⋯\n",
"│ 6.0 │ 2.1239e-6 │ -5.67287 │ 1.68028 │ 0.61638 ⋯\n",
"│ 7.0 │ 3.15775e-10 │ -9.50062 │ 1.67475 │ 0.47672 ⋯\n",
"└───────────┴────────────────┴───────────┴─────────┴────────────────────────────\n",
"\u001b[36m 1 column omitted\u001b[0m\n",
"theoretical values of (α, μ) = (1.618033988749895, 0.5259323034521867)\n"
]
}
],
"source": [
"# g(x) = x^2 - 2 is already defined \n",
"\n",
"a = (1/2)*(1 + sqrt(5));\n",
"mu = ( 2/(2*2*sqrt(2)) )^(a/(1+a));\n",
"\n",
"orderOfConvergence( secant( g, 1. ,2. ), sqrt(2); α=a )\n",
"println(\"theoretical values of (α, μ) = (\", a, \", \" , mu, \")\" )\n"
]
},
{
"cell_type": "markdown",
"id": "22c45f76",
"metadata": {},
"source": [
"Proof. \n",
"\n",
"Again, we may apply Taylor's remainder theorem: there exists $\\eta_n$ between $x_n$ and $x_{n-1}$ such that $f[x_n, x_{n-1}] = f'(\\eta_n)$ and $\\theta_n$ between $x_n$ and $\\xi$ such that $f(x_n) = f'(\\xi)(x_n - \\xi)$. As a result, we have \n",
"\n",
"\\begin{align}\n",
" x_{n+1} - \\xi &= x_n - \\xi - \\frac{f(x_n)}{ f[x_n, x_{n-1}] } \\nonumber\\\\\n",
" %\n",
" &= \\left( 1 - \\frac{f'(\\theta_n)}{ f'(\\eta_n) } \\right) (x_n - \\xi) \n",
"\\end{align}\n",
"\n",
"If $f'(\\xi) > 0$, then there exists a closed interval $I = [\\xi-\\delta, \\xi+\\delta]$ such that \n",
"\n",
"\\begin{align*}\n",
" 0 < \\frac34 f'(\\xi) < f'(x) < \\frac54 f'(\\xi)\n",
"\\end{align*}\n",
"\n",
"for all $x \\in I$. If $x_n, x_{n-1} \\in I$ then $\\eta_n, \\theta_n \\in I$ and so we have $\\frac23 < 1 - \\frac{f'(\\theta_n)}{ f'(\\eta_n) } < \\frac25$ and thus $x_{n+1} \\in I$ with \n",
"\n",
"\\begin{align}\n",
" |x_{n+1} - \\xi| \\leq \\frac23 |x_{n}-\\xi| \\leq \\cdots \\leq \\left(\\frac23\\right)^n |x_1 - \\xi|.\n",
"\\end{align}\n",
"\n",
"That is, the sequence converges at least linearly. Moreover, by (1), we have $|x_{n+1}-\\xi|/|x_n - \\xi| \\to 0$ as $n\\to\\infty$.\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "f1516bc3",
"metadata": {},
"source": [
"Definition (Efficiency index) \n",
"\n",
"The efficiency index is $E := \\alpha^{1/\\theta}$ where $\\theta$ is the number of function evaluations at each step of the iteration. This measures the order of the method *per function evaluation*\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "11962eab",
"metadata": {},
"source": [
"Example. \n",
"\n",
"Newton when $f'(\\xi), f''(\\xi) \\not= 0$, we have $E = 2^{1/2} \\approx 1.41...$ since the order of convergence is $2$ and we need to evaluate $f$ and $f'$ at $x_n$.\n",
"\n",
"Secant when $f'(\\xi) \\not= 0$, we have $E = \\left[\\frac12\\big( 1 + \\sqrt{5} \\big)\\right]^{1/1} \\approx 1.618...$ since the order of convergence is $\\frac12\\big( 1 + \\sqrt{5} \\big)$ and we only require one function evaluation per iteration.\n",
"\n",
"
"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.11.6",
"language": "julia",
"name": "julia-1.11"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}